rm(list = ls())

## Liberaries
library(readxl)
library(tidyverse)
library(coefplot)
library(plyr)
library(effects)
library(dplyr)
library(ggplot2)
library(ggthemes)
library(patchwork)
library(gridExtra)
library(jtools)
library(stargazer)
library(sjPlot)
library(TSstudio)
library(plotly)
library(plyr)
library(MASS)
library(gmodels)
library(ggpubr)


#Load data
Surveydata  <- read_excel (".../Surveydata.xlsx")  ## Copy your pathway here
mydata <- Surveydata 

## drop first observation
mydata <- mydata[-1,]

# Create an ID column
mydata$id <- 1:nrow(mydata)

## Drop incomplete responses
mydata <- mydata[!(mydata$Finished == "False"),]


## Code 1 participants that took the survey after the invasion: 24 Feb, 10 AM, Mountain Standard Time
mydata$postinvasion<-ifelse(mydata$id<276,0,1)
table(mydata$postinvasion)


## Build a character variable Pre vs. Post
mydata$post_ch <-ifelse(mydata$postinvasion =="1","Post","Pre")


## Russian image ordinal coding
mydata$rus_ord <-recode(mydata$Q30, 'Very favorable'=1, 'Somewhat favorable'=2, 
                           'Somewhat unfavorable'=3, 'Very unfavorable'=4)

mydata$rusp_ord <-recode(mydata$Q31, 'Very favorable'=1, 'Somewhat favorable'=2, 
                        'Somewhat unfavorable'=3, 'Very unfavorable'=4)

mydata$rusg_ord <-recode(mydata$Q32, 'Very favorable'=1, 'Somewhat favorable'=2, 
                         'Somewhat unfavorable'=3, 'Very unfavorable'=4)

## Russian image dummy
mydata$rus<-ifelse(mydata$Q30=="Very unfavorable" | mydata$Q30== "Somewhat unfavorable",1,0)

mydata$rusp<-ifelse(mydata$Q31=="Very unfavorable" | mydata$Q31== "Somewhat unfavorable",1,0)

mydata$rusg<-ifelse(mydata$Q32=="Very unfavorable" | mydata$Q32== "Somewhat unfavorable",1,0)



## Russian image dummy: only very unfavorable
mydata$rus_very<-ifelse(mydata$Q30=="Very unfavorable",1,0)
mydata$rusp_very<-ifelse(mydata$Q31=="Very unfavorable",1,0)
mydata$rusg_very<-ifelse(mydata$Q32=="Very unfavorable",1,0)

## Russian people very favorable
mydata$rusp_veryf<-ifelse(mydata$Q31=="Very favorable",1,0)
mydata$rusg_veryf<-ifelse(mydata$Q32=="Very favorable",1,0)
mydata$rus_veryf<-ifelse(mydata$Q30=="Very favorable",1,0)

#gender
mydata$gender<-ifelse(mydata$Q50=="Female",1,0)

#Ideology
mydata$moderate<-ifelse(mydata$Q36_7>3 & mydata$Q36_7<8,1,0)
mydata$conservative<-ifelse(mydata$Q36_7>5,1,0)
mydata$progressive<-ifelse(mydata$Q36_7<6,1,0)

#income
mydata$income.ord<-recode(mydata$Q45, 'Less than $20,000'=1, '$20,000 through $34,999'=1, '$35,000 through $49,999'=2,
'$50,000 through $74,999'=2, '$75,000 through $99,999'=3, '$100,000 or greater'=3)


mydata$income<-ifelse(mydata$Q45 =="$75,000 through $99,999" | mydata$Q45 =="$100,000 or greater",1,0)



# Military family member
mydata$military<-ifelse(mydata$Q22=="Yes",1,0)


#party affiliation
mydata$Democrat<-ifelse(mydata$Q44=="Democrat",1,0)
mydata$Republican<-ifelse(mydata$Q44=="Republican",1,0)
mydata$Independent<-ifelse(mydata$Q44=="Independent",1,0)


#party categorical
mydata$party_categ<-recode(mydata$Q44, 'Democrat'=0, 'Republican'=1, 'Independent'=2)
                   
table(mydata$party_categ)
#Race
mydata$white<-ifelse(mydata$Q39 =="White",1,0)
mydata$black<-ifelse(mydata$Q39=="Black or African American",1,0)



#Political Knowledge
mydata$putin<-ifelse(mydata$Q29=="Vladimir Putin",1,0)
mydata$blinken<-ifelse(mydata$Q34=="Antony Blinken",1,0)
mydata$political_knowledge<- mydata$putin*mydata$blinken


#Religiosity
mydata$religiosity.ord<-recode(mydata$Q21, 'Not at all important'=1, 'Not very important'=2,
                      'Rather important'=3, 'Very important'=4)

mydata$religiosity<-ifelse(mydata$Q21 =="Rather important" | mydata$Q21 =="Rather important",1,0)


##National identity
mydata$national_identity<-ifelse(mydata$Q42 =="Very proud" | mydata$Q42 =="Quite proud",1,0)


## Control and experimental groups
mydata$control<-ifelse(mydata$FL_14_DO =="Control",1,0)
mydata$treatment1<-ifelse(mydata$FL_14_DO =="Stimulus1",1,0)
mydata$treatment2<-ifelse(mydata$FL_14_DO =="Stimulus2",1,0)


## Cponvert to factor
mydata$postinvasion_fac<- as.factor(mydata$postinvasion)
mydata$gender_fac<- as.factor(mydata$gender)
mydata$Democrat_fac<- as.factor(mydata$Democrat)


##Descriptive analysis

#Mean over time
rusm<-aggregate(rus_ord ~ Day, mydata, mean)
ruspm<-aggregate(rusp_ord ~ Day, mydata, mean)
rusgm<-aggregate(rusg_ord ~ Day, mydata, mean)

d1<-merge(x = rusm, y = ruspm, by = c("Day"), all = TRUE)
time_series<-merge(x = d1, y = rusgm, by = c("Day"), all = TRUE)


#country
plot_ly(rusm, x = ~Day, y = ~rus_ord, type = 'scatter', mode = 'lines') %>%
  add_segments(y = 2.7, yend = 3.5, x = as.Date('2022-02-24'), xend = as.Date('2022-02-24'), 
               line = list(dash = "dash"), showlegend=FALSE) %>%
  layout(title = 'Country', plot_bgcolor = "", xaxis = list(title = ''), 
         yaxis = list(title = 'Unfavorability'), legend = list(title=list(text='<b> Species of Iris </b>'))) 



# people
plot_ly(ruspm, x = ~Day, y = ~rusp_ord, type = 'scatter', mode = 'lines') %>%
  add_segments(y = 1.77, yend = 2.1, x = as.Date('2022-02-24'), xend = as.Date('2022-02-24'), 
               line = list(dash = "dash"), showlegend=FALSE) %>%
  layout(title = 'People', plot_bgcolor = "", xaxis = list(title = ''), 
         yaxis = list(title = 'Unfavorability'), legend = list(title=list(text='<b> Species of Iris </b>')))



# Government
plot_ly(rusgm, x = ~Day, y = ~rusg_ord, type = 'scatter',
             mode = 'lines') %>%
  add_segments(y = 3.4, yend = 3.9, x = as.Date('2022-02-24'),  xaxis = list(title = ""),
               xend = as.Date('2022-02-24'), 
               line = list(dash = "dash"), showlegend=FALSE) %>%
  layout(title = 'Government', plot_bgcolor = "", xaxis = list(title = ''), 
         yaxis = list(title = 'Unfavorability'), legend = list(title=list(text='<b> Species of Iris </b>')))






#######Mean scores

#Entire sample
prop.table(table(mydata$rus))

#Pre-intervention sample
prein_sample <- subset(mydata, postinvasion == 0)
prop.table(table(prein_sample$rus))

#Post-intervention sample
postin_sample <- subset(mydata, postinvasion == 1)
prop.table(table(postin_sample$rus))



##Percentage pre vs post

#Pre-intervention sample
prein_sample <- subset(mydata, postinvasion == 0)
prop.table(table(prein_sample$rus))
prop.table(table(prein_sample$rusp))
prop.table(table(prein_sample$rusg))

#Post-intervention sample
postin_sample <- subset(mydata, postinvasion == 1)
prop.table(table(postin_sample$rus))
prop.table(table(postin_sample$rusp))
prop.table(table(postin_sample$rusg))


#Pre-intervention sample: Very unfavorable
prop.table(table(prein_sample$rus_very))
prop.table(table(prein_sample$rusp_very))
prop.table(table(prein_sample$rusg_very))

#Post-intervention sample: Very unfavorable
prop.table(table(postin_sample$rus_very))
prop.table(table(postin_sample$rusp_very))
prop.table(table(postin_sample$rusg_very))




###T-test of means
rus_ttest <- t.test(rus_ord ~ postinvasion, data = mydata, paired = FALSE)
rus_ttest

rusp_ttest <- t.test(rusp_ord ~ postinvasion, data = mydata, paired = FALSE)
rusp_ttest

rusg_ttest <- t.test(rusg_ord ~ postinvasion, data = mydata, paired = FALSE)
rusg_ttest


###T-test of means: Very unfavorable
rus_ttestv <- t.test(rus_very ~ postinvasion, data = mydata, paired = FALSE)
rus_ttestv

rusp_ttestv <- t.test(rusp_very ~ postinvasion, data = mydata, paired = FALSE)
rusp_ttestv

rusg_ttestv <- t.test(rusg_very ~ postinvasion, data = mydata, paired = FALSE)
rusg_ttestv


ggboxplot(mydata, x = "postinvasion", y = "rusg_very", 
          color = "postinvasion", palette = c("#00AFBB", "#E7B800"),
          ylab = "Weight", xlab = "Groups")




### Inferential Analysis

## Attitude toward Russia
mylogit1 <- glm(rus ~ postinvasion_fac + Republican + Democrat  + conservative  +
                 white + black + political_knowledge +
                income + gender + religiosity + military + national_identity 
                , data = mydata, family = binomial())
summary(mylogit1)

ppl1 <- predictorEffect("postinvasion_fac", mylogit1)


## Attitude toward Russian people
mylogit2 <- glm(rusp ~ postinvasion_fac + Republican + Democrat + conservative +
                 white + black + political_knowledge +
                  income + gender + religiosity + military + national_identity, data = mydata, family = binomial())
summary(mylogit2)

ppl2 <- predictorEffect("postinvasion_fac", mylogit2)

## Attitude toward Russian government
mylogit3 <- glm(rusg ~ postinvasion_fac + Republican + Democrat  + conservative +
                   white + black + political_knowledge +
                  income + gender + religiosity + military + national_identity, data = mydata, family = binomial())
summary(mylogit3)

ppl3 <- predictorEffect("postinvasion_fac", mylogit3)




### regression model
pdf(file = "\\reg1.pdf")


plot_summs(mylogit1, mylogit2, mylogit3, coefs = c("Post-invasion" = "postinvasion_fac", "Republican" = "Republican",
                                                  "Democrat" = "Democrat",
                                                   "Conservative" = "conservative", "Moderate" = "moderate",
                                                   "Liberal" = "progressive", "White" = "white", 
                                                   "Black" = "black", "Political knowledge" = "political_knowledge", 
                                                   "Income" = "income", "Gender (female)" = "gender", 
                                                   "Religiosity" = "religiosity", "Military family" = "military", 
                                                  "National pride" = "national_identity"), scale = TRUE, 
           model.names = c("Country", "People", "Government"), legend.title = "Response")

dev.off()



###Predicted probabilities 




#Country
pp1<-sjPlot::plot_model(mylogit1, type = "pred", terms = "postinvasion_fac") +
ggtitle("Country") + theme_set(theme_bw())+
  ylab("") +
  xlab("") + scale_x_continuous(breaks=c(0,1),labels = c("Pre","Post"))

#People
pp2<-sjPlot::plot_model(mylogit2, type = "pred", terms = "postinvasion_fac") +
  ggtitle("People") + theme_set(theme_bw())+
  ylab("") +
  xlab("") + scale_x_continuous(breaks=c(0,1),labels = c("Pre","Post"))

#Government
pp3<-sjPlot::plot_model(mylogit3, type = "pred", terms = "postinvasion_fac") +
  ggtitle("Government") + theme_set(theme_bw())+
  ylab("") +
  xlab("") + scale_x_continuous(breaks=c(0,1),labels = c("Pre","Post"))

pdf(file = "C:\\Users\\pasad\\Dropbox\\Research Projects\\War and image\\predicted1.pdf")
grid.arrange(pp1,pp2,pp3, layout_matrix=rbind(c(1,1,2,2), c(NA, 3, 3, NA)))

dev.off()



##Analysis
mylogit4 <- glm(rus_very ~ postinvasion_fac + Republican + Democrat  + conservative  +
                  white + black + political_knowledge +
                  income + gender + religiosity + military + national_identity
                , data = mydata, family = binomial())


mylogit6 <- glm(rusg_very ~ postinvasion_fac + Republican + Democrat  + conservative  +
                  white + black + political_knowledge +
                  income + gender + religiosity + military + national_identity
                , data = mydata, family = binomial())

### regression model
pdf(file = "\\reg2.pdf")


plot_summs(mylogit4, mylogit6, coefs = c("Post-invasion" = "postinvasion_fac", "Republican" = "Republican",
                                                   "Democrat" = "Democrat",
                                                   "Conservative" = "conservative", "White" = "white", 
                                                   "Black" = "black", "Political knowledge" = "political_knowledge", 
                                                   "Income" = "income", "Gender (female)" = "gender", 
                                                   "Religiosity" = "religiosity", "Military family" = "military", 
                                                   "National pride" = "national_identity"), scale = TRUE, 
           model.names = c("Country", "Government"), legend.title = "Response")

dev.off()


#Predicted probability


#Country
pp4<-sjPlot::plot_model(mylogit4, type = "eff", terms = "postinvasion_fac") +
  ggtitle("Country") + theme_set(theme_bw())+
  ylab("") +
  xlab("") + scale_x_continuous(breaks=c(0,1),labels = c("Pre","Post"))



#Government
pp6<-sjPlot::plot_model(mylogit6, type = "eff", terms = "postinvasion_fac") +
  ggtitle("Government") + theme_set(theme_bw())+
  ylab("") +
  xlab("") + scale_x_continuous(breaks=c(0,1),labels = c("Pre","Post"))


pdf(file = "\\predictedprop.pdf",  width=15, height=5)

ggarrange(
  pp4, pp6, labels = c("", ""),
  common.legend = FALSE, legend = "bottom"
)

dev.off()





## Attitude shift among subgroups towards government

##Create interaction variables


#Female
gendersample <- subset(mydata, gender == 1)

prop.table(table(gendersample$rusg_very))

sjPlot::tab_xtab(var.row = gendersample$post_ch, var.col = gendersample$rusg_very, title = "Table Title", show.row.prc = TRUE)

#Male
mydata$male<-ifelse(mydata$Q50=="Male",1,0)
malesample <- subset(mydata, male == 1)
sjPlot::tab_xtab(var.row = malesample$post_ch, var.col = malesample$rusg_very, title = "Table Title", show.row.prc = TRUE)

#Income 75k and higher
income75sample <- subset(mydata, income == 1)
sjPlot::tab_xtab(var.row = income75sample$post_ch, var.col = income75sample$rusg_very, title = "Table Title", show.row.prc = TRUE)

##Income below 75K
incomebelow75k<- subset(mydata, income == 1)
sjPlot::tab_xtab(var.row = incomebelow75k$post_ch, var.col = incomebelow75k$rusg_very, title = "Table Title", show.row.prc = TRUE)


#Republican
republicansample <- subset(mydata, Republican == 1)
sjPlot::tab_xtab(var.row = republicansample$post_ch, var.col = republicansample$rusg_very, title = "Table Title", show.row.prc = TRUE)


#Democrat
democractsample <- subset(mydata, Democrat == 1)
sjPlot::tab_xtab(var.row = democractsample$post_ch, var.col = democractsample$rusg_very, title = "Table Title", show.row.prc = TRUE)


#Independent
independentsample <- subset(mydata, Independent == 1)
sjPlot::tab_xtab(var.row = independentsample$post_ch, var.col = independentsample$rusg_very, title = "Table Title", show.row.prc = TRUE)


#white
whitesample <- subset(mydata, white == 1)
sjPlot::tab_xtab(var.row = whitesample$post_ch, var.col = whitesample$rusg_very, title = "Table Title", show.row.prc = TRUE)

#Politicak knowledge
pksample <- subset(mydata, political_knowledge == 1)
sjPlot::tab_xtab(var.row = pksample$post_ch, var.col = pksample$rusg_very, title = "Table Title", show.row.prc = TRUE)



#Religiosity
mydata$religiositydummy<-ifelse(mydata$religiosity =="3" | mydata$religiosity =="4" ,1,0)
religiositysample <- subset(mydata, religiosity == 1)
sjPlot::tab_xtab(var.row = religiositysample$post_ch, var.col = religiositysample$rusg_very, title = "Table Title", show.row.prc = TRUE)


#military family
mfsample <- subset(mydata, military == 1)
sjPlot::tab_xtab(var.row = mfsample$post_ch, var.col = mfsample$rusg_very, title = "Table Title", show.row.prc = TRUE)


#national pride
nationalpride <- subset(mydata, national_identity == 1)
sjPlot::tab_xtab(var.row = nationalpride$post_ch, var.col = nationalpride$rusg_very, title = "Table Title", show.row.prc = TRUE)




##Create dataframe for dot plot
Subgroup <- c("Female", "Female", "Male", "Male", "Income > 75k","Income > 75k", "Income < 75k", "Income < 75k",
              "Republican", "Republican", "Democrat", "Democrat", "Independent", "Independent", "White", "White",
              "Politically knowledgeable", "Politically knowledgeable", "Religious", "Religious",  "Military family",
              "Military family", "Proud American", "Proud American")

Group <- c("Pre", "Post", "Pre", "Post", "Pre", "Post", "Pre", "Post", "Pre", "Post", "Pre", "Post", "Pre", "Post",
           "Pre", "Post", "Pre", "Post", "Pre", "Post", "Pre", "Post", "Pre", "Post")

Percent <- c("50.3", "82.3", "62.4", "78.4", "59.6", "80.2", "51.3", "81.1", "39.6", "72.7", "63.5", "82.8",
             "54.1", "81.1", "57.4", "84.5", "58.6", "77.8", "47.6", "80", "65.6", "81", "54.1", "83.3")

df <- data.frame(Subgroup, Group, Percent)

df$Percent <- as.numeric (df$Percent)


## Dot plot
  df %>%
  group_by(Subgroup) %>%
  mutate(Difference = Percent - lead(Percent),
         Position = Percent - (0.5 * Difference)) %>%
  ggplot(aes(x = Percent, y = Subgroup)) +
  geom_vline(xintercept = c(56, 80.6),
             lty = 2, alpha = 0.5) +
  geom_line(aes(group = Subgroup)) +
  geom_point(aes(color = Group), size=3) +
  geom_text(aes(label = abs(Difference),
                y = Subgroup,
                x = Position),
            nudge_y = 0.2) +
  theme_classic(base_size = 12) +
  theme(legend.position = "top") +
  scale_color_brewer(palette = "Paired", direction = -1)



   dplot1<-df %>%
    group_by(Subgroup) %>%
    mutate(Difference = Percent - lead(Percent),
           Position = Percent - (0.5 * Difference)) %>%
    ggplot(aes(x = Percent, y = Subgroup)) +
    geom_vline(xintercept = c(56, 80.6),
               lty = 2, alpha = 0.5) +
    annotate(geom = "text",
             x = c(56),
             y = "",
             label = "Mean: 56") +
    annotate(geom = "text",
             x = c(80.6),
             y = "",
             label = "Mean: 80.6") +
    coord_cartesian(clip = "off") +
    geom_line(aes(group = Subgroup)) +
    geom_point(aes(color = Group), size=3) +
    geom_text(aes(label = abs(Difference),
                  y = Subgroup,
                  x = Position),
              nudge_y = 0.25, size = 3) +
    theme_classic(base_size = 12) +
    theme(legend.position = "top") +
    scale_color_brewer(palette = "Paired", direction = -1) +
    theme(axis.ticks.y = element_line(colour = c("transparent",
                                                 rep("black", 12))))


  
  ##Attitude shift towards country
  
  #gender
  
  #Female
  gendersample <- subset(mydata, gender == 1)
  
  prop.table(table(gendersample$rusg_very))
  
  sjPlot::tab_xtab(var.row = gendersample$post_ch, var.col = gendersample$rus_very, title = "Table Title", show.row.prc = TRUE)
  
  #Male
  mydata$male<-ifelse(mydata$Q50=="Male",1,0)
  malesample <- subset(mydata, male == 1)
  sjPlot::tab_xtab(var.row = malesample$post_ch, var.col = malesample$rus_very, title = "Table Title", show.row.prc = TRUE)
  
  #Income 75k and higher
  income75sample <- subset(mydata, income == 1)
  sjPlot::tab_xtab(var.row = income75sample$post_ch, var.col = income75sample$rus_very, title = "Table Title", show.row.prc = TRUE)
  
  ##Income below 75K
  incomebelow75k<- subset(mydata, income == 1)
  sjPlot::tab_xtab(var.row = incomebelow75k$post_ch, var.col = incomebelow75k$rus_very, title = "Table Title", show.row.prc = TRUE)
  
  
  #Republican
  republicansample <- subset(mydata, Republican == 1)
  sjPlot::tab_xtab(var.row = republicansample$post_ch, var.col = republicansample$rus_very, title = "Table Title", show.row.prc = TRUE)
  
  
  #Democrat
  democractsample <- subset(mydata, Democrat == 1)
  sjPlot::tab_xtab(var.row = democractsample$post_ch, var.col = democractsample$rus_very, title = "Table Title", show.row.prc = TRUE)
  
  
  #Independent
  independentsample <- subset(mydata, Independent == 1)
  sjPlot::tab_xtab(var.row = independentsample$post_ch, var.col = independentsample$rus_very, title = "Table Title", show.row.prc = TRUE)
  
  
  #white
  whitesample <- subset(mydata, white == 1)
  sjPlot::tab_xtab(var.row = whitesample$post_ch, var.col = whitesample$rus_very, title = "Table Title", show.row.prc = TRUE)
  
  #Politicak knowledge
  pksample <- subset(mydata, political_knowledge == 1)
  sjPlot::tab_xtab(var.row = pksample$post_ch, var.col = pksample$rus_very, title = "Table Title", show.row.prc = TRUE)
  
  
  
  #Religiosity
  religiositysample <- subset(mydata, religiosity == 1)
  sjPlot::tab_xtab(var.row = religiositysample$post_ch, var.col = religiositysample$rus_very, title = "Table Title", show.row.prc = TRUE)
  
  
  #military family
  mfsample <- subset(mydata, military == 1)
  sjPlot::tab_xtab(var.row = mfsample$post_ch, var.col = mfsample$rus_very, title = "Table Title", show.row.prc = TRUE)
  
  
  #national pride
  nationalpride <- subset(mydata, national_identity == 1)
  sjPlot::tab_xtab(var.row = nationalpride$post_ch, var.col = nationalpride$rus_very, title = "Table Title", show.row.prc = TRUE)
  
  
  
  
  ##Create dataframe for dot plot
  Subgroup <- c("Female", "Female", "Male", "Male", "Income > 75k","Income > 75k", "Income < 75k", "Income < 75k",
                "Republican", "Republican", "Democrat", "Democrat", "Independent", "Independent", "White", "White",
                "Politically knowledgeable", "Politically knowledgeable", "Religious", "Religious",  "Military family",
                "Military family", "Proud American", "Proud American")
  
  Group <- c("Pre", "Post", "Pre", "Post", "Pre", "Post", "Pre", "Post", "Pre", "Post", "Pre", "Post", "Pre", "Post",
             "Pre", "Post", "Pre", "Post", "Pre", "Post", "Pre", "Post", "Pre", "Post")
  
  Percent <- c("23.5", "35.5", "22.2", "43.1", "24.4", "40.5", "21", "36.8", "24.5", "48.5", "21.2", "39.4",
               "24.7", "34.7", "22", "43.4", "24.7", "39.2", "17.1", "41.1", "32.8", "41.3", "27.7", "45.5")
  
  dfc <- data.frame(Subgroup, Group, Percent)
  
  dfc$Percent <- as.numeric (dfc$Percent)
  

  
  
  
  dplot2 <- dfc %>%
    group_by(Subgroup) %>%
    mutate(Difference = Percent - lead(Percent),
           Position = Percent - (0.5 * Difference)) %>%
    ggplot(aes(x = Percent, y = Subgroup)) +
    geom_vline(xintercept = c(22.9, 38.8),
               lty = 2, alpha = 0.5) +
    annotate(geom = "text",
             x = c(22.9),
             y = "",
             label = "Mean: 22.9") +
    annotate(geom = "text",
             x = c(38.8),
             y = "",
             label = "Mean: 38.8") +
    coord_cartesian(clip = "off") +
    geom_line(aes(group = Subgroup)) +
    geom_point(aes(color = Group), size=3) +
    geom_text(aes(label = abs(Difference),
                  y = Subgroup,
                  x = Position),
              nudge_y = 0.25, size = 3) +
    theme_classic(base_size = 12) +
    theme(legend.position = "top") +
    scale_color_brewer(palette = "Paired", direction = -1) +
    theme(axis.ticks.y = element_line(colour = c("transparent",
                                                 rep("black", 12))))
  
  pdf(file = "\\dplots.pdf", width=9.5, height=6.5)
  
  ggarrange(
    dplot2, dplot1, labels = c("Country", "Government"),
    common.legend = FALSE, legend = "bottom"
  )
  
  dev.off()
  
  
  
    


################################
#############APPENDIX###########
################################


####Balance check

#Gender
mydata$Q50 <- as.factor(mydata$Q50)
ggplot(data=mydata, aes(post_ch))+
  scale_y_continuous(labels=scales::percent) +
  geom_bar(aes(fill=(Q50)),width=0.5, position="fill")+
  labs(x= "",y= "Percent", fill = "Gender")+   scale_colour_wsj("colors6", "") + theme_wsj(color = "gray")+
  coord_flip()



#Income
mydata$Q45 <- as.factor(mydata$Q45)
ggplot(data=mydata, aes(post_ch))+ 
  scale_y_continuous(labels=scales::percent) +
  geom_bar(aes(fill=(Q45)), width=0.5, position="fill")+ 
  labs(x= "",y= "Percent", fill = "Income")+  scale_colour_wsj("colors6", "") + theme_wsj(color = "gray")+
  coord_flip()

##Size   725 * 325

#Party
mydata$Q44 <- as.factor(mydata$Q44)

 ggplot(data=mydata, aes(post_ch))+
  scale_y_continuous(labels=scales::percent) +
  geom_bar(aes(fill=(Q44)), width=0.5, position="fill")+
  labs(x= "",y= "Percent", fill = "Party")+  theme_tufte(ticks = FALSE)+  theme_wsj(color = "gray")+
  coord_flip()

#Race
mydata$Q38...54 <- as.factor(mydata$Q38...54)

 ggplot(data=mydata, aes(post_ch))+
  scale_y_continuous(labels=scales::percent) +
  geom_bar(aes(fill=(Q38...54)),width=0.5, position="fill")+
  labs(x= "",y= "Percent", fill = "Race")+  theme_tufte(ticks = FALSE)+ theme_wsj(color = "gray")+
  coord_flip()

#Political knowledge
mydata$Q29 <- as.factor(mydata$Q29)

ggplot(data=mydata, aes(post_ch))+
  scale_y_continuous(labels=scales::percent) +
  geom_bar(aes(fill=(Q29)), width=0.5, position="fill")+
  labs(x= "",y= "Percent", fill = "President")+  theme_tufte(ticks = FALSE)+ theme_wsj(color = "gray")+
  coord_flip()


mydata$Q34 <- as.factor(mydata$Q34)

ggplot(data=mydata, aes(post_ch))+
  scale_y_continuous(labels=scales::percent) +
  geom_bar(aes(fill=(Q34)),width=0.5, position="fill")+
  labs(x= "",y= "Percent", fill = "Secretary ")+  theme_tufte(ticks = FALSE)+ theme_wsj(color = "gray")+
  coord_flip()

#Religiosity
mydata$Q21 <- as.factor(mydata$Q21)

 ggplot(data=mydata, aes(post_ch))+
  scale_y_continuous(labels=scales::percent) +
  geom_bar(aes(fill=(Q21)),width=0.5, position="fill")+
  labs(x= "",y= "Percent", fill = "Religiosity") + theme_tufte(ticks = FALSE)+ theme_wsj(color = "gray")+
  coord_flip()


#Military
mydata$Q22 <- as.factor(mydata$Q22)

ggplot(data=mydata, aes(post_ch))+
  scale_y_continuous(labels=scales::percent) +
  geom_bar(aes(fill=(Q22)),width=0.5, position="fill")+
  labs(x= "",y= "Percent", fill = "Military family member") + theme_tufte(ticks = FALSE)+ theme_wsj(color = "gray")+
  coord_flip()


#ideology
mydata$Q36_7 <- as.factor(mydata$Q36_7)

ggplot(data=mydata, aes(post_ch))+
  scale_y_continuous(labels=scales::percent) +
  geom_bar(aes(fill=(Q36_7)),width=0.5, position="fill")+ 
  labs(x= "",y= "Percent", fill = "Ideology")+  theme_wsj(color = "gray")+
  coord_flip()

##National pride
mydata$Q42 <- as.factor(mydata$Q42)
ggplot(data=mydata, aes(post_ch))+ 
  scale_y_continuous(labels=scales::percent) +
  geom_bar(aes(fill=(Q42)), width=0.5, position="fill")+
  labs(x= "",y= "Percent", fill = "National pride")+  scale_colour_wsj("colors6", "") + theme_wsj(color = "gray")+
  coord_flip()

#####T-test of difference in means#########


#Gender
group_by(mydata, postinvasion) %>%
  summarise(
    count = n(),
    mean = mean(gender, na.rm = TRUE),
    sd = sd(gender, na.rm = TRUE)
  )

gendertt <- t.test(gender ~ postinvasion, data = mydata, paired = FALSE)
gendertt


#Income
incomett <- t.test(income ~ postinvasion, data = mydata, paired = FALSE)
incomett

incomemw<-wilcox.test(income ~ postinvasion, data=mydata, na.rm=TRUE, paired=FALSE, exact=FALSE, conf.int=TRUE)
print(incomemw)

#Race
racett <- t.test(white ~ postinvasion, data = mydata, paired = FALSE)
racett

racettb <- t.test(black ~ postinvasion, data = mydata, paired = FALSE)
racettb

#Party

reptt <- t.test(Republican ~ postinvasion, data = mydata, paired = FALSE)
reptt

demtt <- t.test(Democrat ~ postinvasion, data = mydata, paired = FALSE)
demtt

#Ideology
ideott <- t.test(as.numeric(Q36_7) ~ postinvasion, data = mydata, paired = FALSE)
ideott

ideomw<-wilcox.test(as.numeric(Q36_7) ~ postinvasion, data=mydata, na.rm=TRUE, paired=FALSE, exact=FALSE, conf.int=TRUE)
print(ideomw)


#Political knowledge
pktt <- t.test(political_knowledge ~ postinvasion, data = mydata, paired = FALSE)
pktt

#Religiosity
relitt <- t.test(religiosity ~ postinvasion, data = mydata, paired = FALSE)
relitt

#Military
miltt <- t.test(military ~ postinvasion, data = mydata, paired = FALSE)
miltt

#National pride
nationaltt <- t.test(national_identity ~ postinvasion, data = mydata, paired = FALSE)
nationaltt 




##Ordered logit
myordered1<-polr (as.factor(rus_ord) ~ postinvasion_fac + Republican + Democrat  + conservative  +
      white + black + political_knowledge +
      income + gender + religiosity + military + national_identity 
    , data = mydata, Hess = TRUE, method = c("logistic"))
summary(myordered1)



myordered2<-polr (as.factor(rusp_ord) ~ postinvasion_fac + Republican + Democrat  + conservative  +
                    white + black + political_knowledge +
                    income + gender + religiosity + military + national_identity 
                  , data = mydata, Hess = TRUE, method = c("logistic"))
summary(myordered2)



myordered3<-polr (as.factor(rusg_ord) ~ postinvasion_fac + Republican + Democrat  + conservative  +
                    white + black + political_knowledge +
                    income + gender + religiosity + military + national_identity 
                  , data = mydata, Hess = TRUE, method = c("logistic"))
summary(myordered3)


pdf(file = "\\orderedlogit.pdf")


plot_summs(myordered1, myordered2, myordered3, coefs = c("Post-invasion" = "postinvasion_fac1", "Republican" = "Republican",
                                                         "Democrat" = "Democrat",
                                                         "Conservative" = "conservative", "White" = "white", 
                                                         "Black" = "black", "Political knowledge" = "political_knowledge", 
                                                         "Income" = "income", "Gender (female)" = "gender", 
                                                         "Religiosity" = "religiosity", "Military family" = "military", 
                                                         "National pride" = "national_identity"), scale = TRUE,
           model.names = c("Country", "People", "Government"), legend.title = "Response")

dev.off()


##PEW and GALLUP
#Data
Favorable <- c("21", "22", "33", "21", "27") ##Percentage
Unfavorable<- c("71", "77", "67", "79", "73")  ##Percentage 
sample <- c("PEW", "Gallup", "Pre", "Post", "Entire sample")
comparison <- data.frame(Favorable , Unfavorable, sample)

##code

comparison %>% 
  pivot_longer(
    cols = -sample,
    names_to = "Category",
    values_to = "percentage", 
    values_transform = list(percentage = as.integer)
  ) %>% 
  ggplot(aes(x = sample, y=percentage, fill=Category))+
  geom_col(position = position_fill()) +
  scale_y_continuous(labels = scales::percent) +
  geom_text(aes(label = percentage),
            position = position_fill(vjust = .3))+ 
  theme_light()

##Regression analysis: Interactions

#Create interaction variables
mydata$gender.int <- mydata$gender*mydata$postinvasion
mydata$republican.int <- mydata$Republican*mydata$postinvasion
mydata$democrat.int <- mydata$Democrat*mydata$postinvasion
mydata$conservative.int <- mydata$conservative*mydata$postinvasion
mydata$white.int <- mydata$white*mydata$postinvasion
mydata$black.int <- mydata$black*mydata$postinvasion
mydata$income.int <- mydata$income*mydata$postinvasion
mydata$relig.int <- mydata$religiosity*mydata$postinvasion
mydata$national.int <- mydata$national_identity*mydata$postinvasion
mydata$pk.int <- mydata$political_knowledge*mydata$postinvasion


##regression
mydata$russiangov_fav<-ifelse(mydata$Q32=="Very favorable" | mydata$Q32== "Somewhat favorable",1,0)
table(mydata$russiangov_fav)


